Synchronization of a Josephson junction array in terms of global variables 
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We consider an array of Josephson junctions with a common LCR-load. Application of the 
Watanabe-Strogatz approach [Physica D, v. 74, p. 197 (1994)] allows us to formulate the dynamics 
of the array via the global variables only. For identical junctions this is a finite set of equations, 
analysis of which reveals the regions of bistability of the synchronous and asynchronous states. For 
disordered arrays with distributed parameters of the junctions, the problem is formulated as an 
integro-differential equation for the global variables, here stability of the asynchronous states and 
the properties of the transition synchrony- asynchrony are established numerically. 

PACS numbers: 05.45.Xt,74.81.Fa 



I. INTRODUCTION 

Synchronization in populations of coupled oscillators 
is a general phenomenon observed in many physical sys- 
tems, see recent experimental studies of optomechanical, 
micromechanical, electronic, mechanical, chemical oscil- 
lators jl| . Synchronization effects are also ubiquitous in 
biology and social sciences. One of the basic examples of 
oscillating physical systems that being coupled synchro- 
nize, are Josephson junctions [2|. In theoretical stud- 
ies of the Josephson junction arrays one typically either 
performs direct numerical simulation of the microscopic 
equations (see, e.g., 0]) or reduces the problem to the 
standard Kuramoto-type model [l-la|. 

Quite remarkable in this respect is the paper [7| , where 
a careful comparison of the microscopic modeling and 
the reduced Kuramoto-type model has been performed. 
The authors demonstrated that a hysteretic transition 
to synchrony in an array of Josephson junctions can be 
explained by a Kuramoto-type modeling (where usually 
the transition is not hysteretic), if in its derivation one 
self-consistently accounts for changes of the oscillator pa- 
rameters. 

Our aim in this paper is to shed light on the hysteretic 
transitions to synchrony in Josephson arrays by studying 
the equations for global variables. In this approach, that 
is based on the seminal papers by Watanabe and Stro- 
gatz (WS) [a, 0; it is possible to formulate exact low- 
dimensional equations for the array, without using ap- 
proximate reduction to the Kuramoto model. The paper 
is organized as follows. First, we formulate the equations 
for the array of identical junctions via the global vari- 
ables. Analysis of these equations shows regions of bista- 
bility asynchrony-synchrony, and the hysteretic transi- 
tions. Then we proceed to non-identical junctions, where 
the equations are of more complex form. Here we analyze 
stability of asynchronous states, and show numerically 
that the transition to synchrony is also hysteretic. 



II. IDENTICAL JUNCTIONS 

A. Formulation in terms of global variables 

We start with formulating the system of equations for 
the Josephson junction series array with a LCR load. 
Our setup is the same as in refs. [j-Q, the equatons for 
the junction phases (pi and the load capacitor charge Q 
read 



h d(pi 
2erlU 
d^Q 



Ic sin ipi = I 



dQ 
dt 



„dQ ,Q___^sr^ ^^ 



(1) 



dt C 2e 



Here N is the number of junctions, described by a resis- 
tive model with critical current I^ and resistance r, while 
L, C, R are parameters of the LCR-load. It is convenient 
to introduce dimensionless variables according to 

LOc — 2erlc/fi^ t* = ujct, 
Q*=u,L*Q/h, r=I/Ic, (2) 

R* = R/rN, L* ^ UJcL/rN, C* = NtJavC , 

and to rewrite the system ([T]) in a dimensionless form 
(droping asterixes for simplicity) 



ip, ^ I - eQ ~ simp,, 

1 ^ 
7Q + ujqQ = / - — ^ sin (/?,, 



(3) 



where e = 1/L*, 7 = (i?* + 1)/L*, and ujq = 1/VL*C*. 

The global coupling can be represented through the 
complex mean field (Kuramoto order parameter) 
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and the equations for the junction phases can be written 
as 



•^i 



I ~eQ + l-m{e-''f^). 



(5) 



This form of the phase equation allows us to use the 
Watanabe-Strogatz ansatz [^, ^, apphcable to general 
systems of phase equations driven by a common force 
and having form 



v3, = /(i)+Im(G(i)e-''^') 



(6) 



with arbitrary real f{t) and complex G{t) (in our case 
f = I — eQ, G = 1). We use the formulation of the 
Watanabe-Strogatz theory presented in Ref. [101 • The 
ensemble is characterized by three global time-dependent 
WS variables p, $, ^ and N constants of motion ipi (of 
which only iV — 3 are independent) which are related to 
the phases ipi as 



, p + exp(i(V>, - ^)) 
pexp(i(V', -*)) + ! 



(7) 



with additional conditions ^^cosV'i = ^^smipi = 
J2iCOs2ijJi = 0. The equations for the global WS vari- 
ables read ^^ 
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To close the system we need to add the equation for Q, 
where the imaginary part of the order parameter Z en- 
ters, so Z should be represented through the Watanabe- 
Strogatz variables. In general, the expression for Z is 
rather complex (cf. [iOillll) but in the case of a uniform 
distribution of the constants tpi, the order parameter is 
just Z = yoe'*. This important case, where WS global 
variables p, $ have a clear physical meaning as the com- 
ponents of the Kuramoto order parameter, will be treated 
below. Additionally, we notice that the variable ^ does 
not enter other equations, so we obtain a closed system 
of equations that describes the array 



Q + 7Q + LJoQ 



i{I-eQ)Z 
I -Im{Z). 



Z^ 
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conditions, we discuss its relevance below. Second, the 
order parameter Z does not vanish in the case of full 
asynchrony of junctions: for noncoupled junctions with 
e = we get a steady state Zq — i(J — \/P — 1). This 
non-vanishing value appears because free junctions ro- 
tate non-uniformly and the "natural" distribution of the 
phases in the asynchronous state is not uniform. 

We start the analysis of © with finding its steady 
states. Because at such a state Q — 0, the coupling van- 
ishes and the steady state describing the asynchronous 
regime with Zq = i(/ — V/^ — 1), Qo ~ uJq'^VP — 1 is 
the only stationary solution. Stability of this solution is 
determined by the fourth-order characteristic equation 



+ [(7 - e)(/2 _ 1) + el^P - 1 ]A + c.2(/2 - 1) = 0. 

(10) 
The stability border can be easily found by assuming 
A = icj: 



ujI = {P - 1) + -Vp - 1(1 - Vp - 1). (11) 

7 

The fully synchronous solution of ([9]) corresponds to 
the case \Z\ — 1, so that only the phase $ changes, ac- 
cording to the system 



Q + 7(5 + w^Q = /-sin$ 
$ = / - eQ - sin $ . 



(12) 



We have found the limit cycle in Eq. (|12p numerically and 
determined its stability by finding the largest multiplier. 
Together with expression (|TT|) this allows us to find the 
domains of stability of the asynchronous and synchronous 
states, together with the region of bistability of these 
regimes, see Fig.[T] 

In Fig. [2] we give another illustration of the bistabil- 
ity, presenting the dependence of Zq on parameter /, to- 
gether with the value \Z\ = 1 in the synchronous case. 
Here we also show what happens if our basic assumption 
at derivation of eqs. ([5]), namely of a uniform distribu- 
tion of constants 4'ii is not satisfied. We have simulated 
an ensemble of 100 junctions, preparing the initial con- 
ditions with a nonuniform distribution of constants ipi 
as described in ref. ^101, appendix C. Instead of leading 
to a stable state Zq, the desynchronous population now 
shows an oscillating variable Z(t), minima and maxima 
of which are marked with squares. In the synchronous 
regime, \Z\ = 1 as before, and the information on the 
constants Vi gets lost as synchrony establishes. 



B. Bistability and hysteretic transitions 

Analysis of system (0) is our goal in the rest of this 
section. Before proceeding, some remarks are in order. 
First, in the derivation of (jH]) no approximation except 
for an assumption of a uniform distribution of constants 
■04, has been made. The latter is a restriction on initial 



III. NONIDENTICAL JUNCTIONS 

A. Formulation of the model 

There are two parameters of individual junctions that 
can differ: the critical current Ic and the resistance r 
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FIG. 1. (Color online) Domains of stability of synchronous 
(above lower dashed line) and asynchronous (below upper 
solid line) states on the plane of parameters {luq,Q^), where 
Q, — y/P — 1 is the natural frequency of the junctions. Here 
e = 0.5, and 7 = 1.0 (a), 1.7 (b), 2.7 (c). 
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FIG. 2. Dependence of the order parameter \Z\ on the current 
I for 100 junctions. Line: uniform distribution of constants 
tpi, squares: nonuniform distributions. 



the junctions read 



t/5fcj = (1 + T]k)[I - eQ - (1 + £,k)sm(pki 
M p 



-iQ + ulQ^I-- ^(1 + r?fc)(l + 6) Y. ^i" 



^ki 



k=l 



(13) 

To each group the Watanabe-Strogatz ansatz as de- 
scribed in the previous section can be applied, and as 
a result instead of the identical array equations ^ we 
obtain a system 



Q + -iQ + ujlQ = I-{{l + 77fe)(l + 6)Ini(Zfe)), 
Zk = {l + Vk) (i{I - eQ)Zk + (1 + S.k)^—7^ 



(14) 



where average () is taken over all groups. Starting from 
(|14p one can easily take a thermodynamic limit of an in- 
finite number of groups M — >■ 00, in this limit Z^ — >■ 
Z{r],^). Then ([T4|) reduces to an integro-differential 
equation that includes the distribution function W(r],^) 
of disorder parameters ^,77 (cf. [lO.]): 



(cf. fs, i>|). In order to be able to apply the Watanabe- 
Strogatz approach as above, we will assume that they are 
organized in groups, each of the size P, and the parame- 
ters of all junctions in a group are identical: the critical 
current is /c(l + Cfc) and the resistance is r{l + rik), where 
index k = 1, . . . , Af counts the groups. The total number 
of junctions is A^ = MP. In this setup the equations for 



Q + -/Q + cjIQ = I- 
d77 dC W^(77, (1 + r/) (1 + e)Ini(Z(r7, 0) , 

1-Z2 



ziv,0^a + v)[Ki-(Q)z + ii + 0- 



(15) 




FIG. 3. (Color online) Real part of the maximum eigenvalue 
A as a function of the dimensionless current / for the different 
values (numbers on the panels) of fi (panel (a)) and ( (panel 
(b)). 



B. Asynchronous state and its stability 

The asynchronous state is the steady state of the sys- 
tem (113 : 



^o(?7,0 =i- 



1+^ 



(16) 
where we assume {£) = (rj) = 0. Remarkably, the dis- 
order in the junction resistances (parameter 77) does not 
influence the value Zg, only the disorder in critical cur- 
rents (parameter ^). However, the stability of this asyn- 
chronous state depends on distributions of 77 and ^. We 
consider two cases, with a disorder in one parameter only, 
(i) Disorder in resistances. Here we assume that 
^{VtO = HO^i^iv) where W^ is a uniform distribu- 
tion in the interval (— /i,^). To study the perturbations 
in the integral equation (fT5l) at the steady solution ([T6| , 
we discretized the integral using 500 nodes and found the 
eigenvalues of the resulting matrix. The results for the 
maximal eigenvalue are shown in Fig. |3^. One can see 
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FIG. 4. (Color online) Panels (a),(d): Dependence of the av- 
eraged order parameter \z\ on current I, fj, — 0.01, C — and 
H = 0, ^ = 0.05 respectively. Three lines show the maximal 
(upper dashed line), the average (solid line), and the minimal 
(lower dashed line) value of variations of \z\ in time, in the 
asynchronous states these lines coincide. Panels (c), (d), (e) 
and (f ) show enlargements of the regions near the synchrony- 
asynchrony transitions, to demonstrate the hysteresis. 



that, with increasing the external current /, the asyn- 
chronous state loses stability almost at the same critical 
value as for identical junctions (expression ((TT|) ). but for 
large values of / the stability is restored. The region of 
instability decreases for larger disorder fj,. 

(ii) Disorder in critical currents. Here we assume that 
^iVi = '5(^)T^c(C)' where ( is the width of the uniform 
distribution. With the same procedure as in case (i) we 
found the stability eigenvalues that are shown in Fig.[3]3. 
Qualitatively, the pictures look similar: both disorders 
result in a finite (in therms of the external current /) 
region of instability of the asynchronous state. 

Both calculations presented in Fig. [3] show, that the 
main effect of disorder in arrays is in the establishing 
of stability of the asynchronous state for large values of 
current /, while only in some range (which decreases with 
disorder) the asynchrony is unstable. We illustrate the 
appearing synchrony patterns in disordered arrays in the 
next subsection. 



C. Numerical simulations 

Dynamics of the nonhomogeneous arrays of Josephson 
junctions is illustrated in Figs. H) As above, we consider 
not a general situation where both the critical current 
and the resistance are spread, but cases where one of 
these parameters has a distribution. In numerical simula- 
tions we use the discrete representation (J14p . In order to 



avoid spurious non-smooth solutions, an additional very 
small viscous term ~ (^fe+i + Zk_i — 2Z^.) was added to 
the equation for Zk that ensures numerical stabilization 
of the integro-differential equation. 

To characterize synchrony we calculated the average 
over the array order parameter z — M~^ ^^ Zk and plot 
it vs. parameter / in Fig. 21 In the asynchronous state 
this parameter attains the fixed point (cf. Eq. (J16l) ). 
while in the synchronous state it oscillates arround some 
mean value (because of disorder the synchrony is not 
complete, so \z\ < 1). Remarkably, also in the case of 
disorder, the transition to synchrony demonstrates hys- 
teresis both for small and large values of /, as can be seen 
on panels (b),(c),(e), and (f) of Fig. S) 



IV. CONCLUSION 

In this paper we applied the approach by Watan- 
abe and Strogatz to the description of the synchro- 
nization transition in an array of Josephson junctions 
with an LCR load. For identical junctions a closed 
low-dimensional system of equations for global variables 
(the Watanabe-Strogatz variables for the junctions and 
two variables describing the load) demonstrates a re- 
gion of bistability at the transition from asynchrony to 



full synchrony, so that this transition shows hysteresis. 
This confirms previous results based on the approximate 
self-consistent reduction to the Kuramoto model Q- 
For nonidentical junction the method yields an integro- 
differential system, as each group of junctions having cer- 
tain parameters is described by the WS variables. Here, 
with the growth of the variability of parameters, the re- 
gion of synchronization shrinks. Transition to synchrony 
in this case is also hysteretic. 

Validity of the WS approach to the Josephson junction 
array is based on the fact, that for standard junctions 
the dependence of the superconducting current on the 
phase is a simple sine function. Therefore, the theory is 
also valid for so-called 7r-junctions [i3| , where the current 
has an opposite direction but nevertheless is propotional 
to sm{ip). However, for recently constructed so-called ip- 
j unctions [1J|, where the phase dependence of the current 
contains the second harmonics, the WS approach is not 
applicable, and synchronization of such junctions remains 
a challenging problem. 
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